Systematic Conservation Planning is a field of research that uses GIS and Python for spatial analysis, with the purpose of strategically identifying areas that meet a given conservation goal.
Scientists bring together collected data from a variety of biological and ecological sources and execute an iterative workflow to define areas that successfully meet the required criteria. These data sources can include:
Results of the analysis are shared with key stakeholders and decision makers, enabling them to efficiently direct limited resources towards solutions with the lowest cost and greatest chance for long term ecological success.
This is a critical time for conservation, as many ecosystems across the globe are at risk of being lost due to deforestation, construction, large-scale agriculture, human encroachment, climate change etc. There are multiple agencies and organizations working on both a regional and global level to fast-track conservation efforts.
In 2021, U.S. President Biden launched the ‘America the Beautiful’ initiative, with the stated goal of conserving at least 30 percent of U.S. lands and waters by 2030. Globally, the International Union for Conservation of Nature (IUCN) has been working since 2004 on a standard for identifying Key Biodiversity Areas (KBAs), which are defined as areas essential to supporting the persistence of global biodiversity at either a species or ecosystem level. Once KBAs are established and the information is shared with stakeholders, that area becomes prioritized for conservation planning.
The IUCN is aware their KBA standard will evolve over time, and actively seeks input from users who work with the current criteria and may have suggestions for improvement. Our project sponsor Pat Comer is a member of the IUCN Commission on Ecosystem Management, and is interested in whether the current thresholds adequately allow for KBAs to be defined for all ecosystems.
Currently, the standard is set so that a KBA can be identified at a site that comprises 10% of the global extent of an ecosystem defined as ‘Vulnerable’ by the IUCN Red List of Ecosystems (RLE), or at 5% when the ecosystem is at the higher risk ‘Critical’ or ‘Endangered’ level.
We would like to investigate if these current thresholds may be set too high, especially for ecosystems with linear or fragmented spatial characteristics. In these instances, it can be difficult to locate independent areas large enough to qualify for conservation using the current standard, or that are in close enough proximity to other areas to allow them to be managed together as a single conservation unit.
The question we are researching is whether lower thresholds would be more successful in identifying KBAs for ecosystems with such spatial characteristics. If so, that may be a useful revision for the IUCN to make to their standard, as ultimately the goal is to identify KBAs for all threatened ecosystems.
Currently our study area includes continental US. For the trial run we used nine ecosystems with diverse spatial characteristics:
# Import libraries
import os
from glob import glob
import matplotlib.pyplot as plt
import numpy as np
import geopandas as gpd
import xarray as xr
import rioxarray as rxr
from rasterio.plot import plotting_extent
import earthpy as et
# Set working directory
os.chdir(os.path.join(et.io.HOME,
'earth-analytics',
'data'))
# Study Area
study_area = gpd.read_file(os.path.join('natureserve',
'study_area',
'us_bound_albers.shp'))
# Set paths to hexes
path_hex = os.path.join('natureserve', 'shp_hex')
ecosytem_hex_paths = sorted(glob(os.path.join(path_hex, '*.shp')))
# Plot the Study Area
fig, ax = plt.subplots(figsize=(20, 12))
ax.set_title('Study Area and Nine Ecosystems Selected for the Trail Run',
fontsize=20)
for afile in ecosytem_hex_paths:
shp = gpd.read_file(afile)
shp.plot(ax=ax, color='blue', zorder=2)
study_area.plot(ax=ax, facecolor='lavender', zorder=1)
ax.axis('off')
(-2586827.152087886, 2488867.850089356, 115760.23776860186, 3310957.9078612914)
# Set paths to raster data
path_raster = os.path.join('natureserve', 'raster')
# In case if all rasters in a folder need to be read/plotted
ecosytem_rasters = sorted(glob(os.path.join(path_raster, '*.tif')))
# Read individual raster
tallgrass_r = rxr.open_rasterio(os.path.join(path_raster,
'tallgrass_r.tif')).squeeze()
# Create mask to visalize nodata raster cells
tallgrass_r_ma = (
np.ma.masked_where(
tallgrass_r != tallgrass_r.rio.nodata, tallgrass_r, copy=True))
# Visualize individual ecosystem raster
fig, ax = plt.subplots(figsize=(20, 12))
ax.set_title('Central Tallgrass Prarie, Raster', fontsize=20)
ax.set(xlabel='Meters',
ylabel='Meters')
plot = ax.imshow(tallgrass_r_ma, cmap='summer')
plt.show()
# Read individual raster
prarie_r = rxr.open_rasterio(os.path.join(path_raster,
'prarie_r.tif')).squeeze()
# Create mask to visalize nodata raster cells
prarie_r_ma = (
np.ma.masked_where(
prarie_r != prarie_r.rio.nodata, prarie_r, copy=True))
# Visualize individual ecosystem raster
fig, ax = plt.subplots(figsize=(20, 12))
ax.set_title('Columbia Basin Palouse Prairie, Raster', fontsize=20)
plot = ax.imshow(prarie_r_ma)
ax.set(xlabel='Meters',
ylabel='Meters')
plt.show()
# Read individual planning unit shape file
tallgrass_v = gpd.read_file(os.path.join('natureserve',
'shp_hex',
'tallgrass_hex.shp'))
# Visualize individual planning unit shape file
fig, ax = plt.subplots(figsize=(20, 12))
ax.set_title('Central Tallgrass Prarie, Planning Units', fontsize=20)
ax.set(xlabel='Meters',
ylabel='Meters')
tallgrass_v.plot(ax=ax, color='none', edgecolor='green')
<AxesSubplot:title={'center':'Central Tallgrass Prarie, Planning Units'}, xlabel='Meters', ylabel='Meters'>
prarie_v = gpd.read_file(os.path.join('natureserve',
'shp_hex',
'prarie_hex.shp'))
fig, ax = plt.subplots(figsize=(20, 12))
ax.set_title('Columbia Basin Palouse Prairie, Planning Units', fontsize=20)
ax.set(xlabel='Meters',
ylabel='Meters')
prarie_v.plot(ax=ax, color='none', edgecolor='purple')
<AxesSubplot:title={'center':'Columbia Basin Palouse Prairie, Planning Units'}, xlabel='Meters', ylabel='Meters'>
# Set plotting extent
tallgrass_extent = plotting_extent(tallgrass_r, tallgrass_r.rio.transform())
# Visualize overlay of individual ecosystem raster and planning units
fig, ax = plt.subplots(figsize=(22, 14))
ax.set_title('Central Tallgrass Prarie, Overlaid With Planning Units',
fontsize=20)
plot = ax.imshow(tallgrass_r_ma, cmap='summer', extent=tallgrass_extent)
tallgrass_v.plot(ax=ax, color='none', edgecolor='yellow')
ax.axis('off')
plt.show()
# Set plotting extent
prarie_extent = plotting_extent(prarie_r, prarie_r.rio.transform())
prarie_extent
# Visualize overlay of individual ecosystem raster and planning units
fig, ax = plt.subplots(figsize=(22, 14))
ax.set_title('Columbia Basin Palouse Prairie, Overlaid With Planning Units',
fontsize=20)
plot = ax.imshow(prarie_r_ma, extent=prarie_extent)
prarie_v.plot(ax=ax, color='none', edgecolor='cyan')
ax.axis('off')
plt.show()
tallgrass_marxan = gpd.read_file(os.path.join('natureserve',
'marxan_prep',
'tallgrass',
'pulayerfeatures.shp'))
fig, ax = plt.subplots(figsize=(22, 14))
ax.set_title('Prarie Surface Area (sq. m.) which Falls within \
Each Planning Unit', fontsize=20)
tallgrass_marxan.plot(column='VALUE_7421', ax=ax, cmap='summer', legend=True)
ax.axis('off')
(-405727.24308490753, 861196.5245847702, 1705617.5982001543, 2576009.613000083)
# Visualize results
prarie_marxan = gpd.read_file(os.path.join('natureserve',
'marxan_prep',
'prarie',
'pulayerfeatures.shp'))
prarie_marxan
fig, ax = plt.subplots(figsize=(22, 14))
ax.set_title('Prarie Surface Area (sq. m.) which Falls within \
Each Planning Unit', fontsize=20)
prarie_marxan.plot(column='VALUE_7142', ax=ax, legend=True)
ax.axis('off')
(-1810739.8217100264, -1515885.191089952, 2681195.5020201444, 2912384.8005802394)
We are aiming to conduct a sensitivity analysis of the current thresholds set by the IUCN for KBA identification. This will allow us to see if a lower threshold will be more effective in identifying areas to prioritize for conservation.
If we find lower thresholds are helpful in identifying KBAs for ecosystems with linear or fragmented spatial characteristics, this could be presented to the IUCN for inclusion in a future revision to the KBA standard which would lead to our findings being used globally.
Systematic conservation planning relies on a complex annealing algorithm function that accounts for a number of variables. This algorithm can be applied to a variety of conditions, depending on the input factors. The question of how to determine which input factors return the best outcomes can be capably answered using an iterative Python workflow to help manage the complexities and analyze the outcomes. An additional benefit of our Python workflow is that it will enable interested parties to review our results in an objective, repeatable and transparent way, or even use their own source data for a similar sensitivity analysis.
We are left with several questions at the end of this initial project phase:
Most obviously, how can we improve the workflow? Any open-source scientific workflow has room for improvement.
Can the entire workflow be completed in Python, or is a GIS software necessary? If it can, how? What tools are necessary?
How will the spatial characteristics of the ecosystems be analyzed affect their calculated KBA status?